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Abstract. In the context of the present and future Cosmic Microwave Background (CMB) experiments, going 
beyond the information provided by the power spectrum has become necessary in order to tightly constrain the 
cosmological model. The non-Gaussian signatures in the CMB represent a very promising tool to probe the early 
universe and the structure formation epoch. We present the results of a comparison between two families of non- 
Gaussian estimators: The first act on the wavelet space (skewness and excess kurtosis of the wavelet coefficients) 
and the second group on the Fourier space (bi- and trispectrum) . We compare the relative sensitivities of these 
estimators by applying them to three different data sets meant to reproduce the majority of possible non-Gaussian 
contributions to the CMB. We find that the skewness in the wavelet space is slightly more sensitive than the 
bispectrum. For the four point estimators, we find that the excess kurtosis of the wavelet coefficients has very 
similar capabilities than the diagonal trispectrum while a near-diagonal trispectrum seems to be less sensitive to 
non-Gaussian signatures. 
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1. Introduction 

Over the last few years, the amount of available data 
on the Cosmic Microwave Background (CMB) has in- 
creased dramatically. Recent measurements of the CMB 
power spectrum (BOOMERANG, MAXIMA, DASI, CBI, 
ARCHEOPS, VSA, ACBAR...) have started to shed light 
on the origin and large scale structure of our universe. 
In particular, the determination of the density parame- 
ter f2 ~ 1 has been obtained from the position of the 
first acoustic peak. Together with results from other ex- 
periments, which probe e.g. the distribution of galaxies 
and the distances to typc-Ia supcrnovae, a coherent im- 
age of todays universe has started to emerge (see e.g. 
Wang et al.(2002)). Testing the detailed statistical nature 
of the CMB anisotropics has become not only one of the 
major goals of CMB cosmology but it is now also within 
our reach. The MAP satellite (and also the future Planck 
mission) will provide us with an all sky survey of the CMB 
with good angular resolution that can be used to search 
for non-Gaussian signatures. 

A Gaussian random field is completely determined 
by its power spectrum. But many models of the 
early universe, like inflation (e.g. Bartolo et al.(2002), 
Bernardeau & Uzan(2002) and references therein), su- 
per strings or topological defects, predict non-Gaussian 
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contributions to the initial fluctuations Luo(1994b), 
Jaffe(1994), Gangui et al.(1994), . Furthermore, any non- 
linearity in their evolution introduces additional non- 
Gaussian signatures. An extreme example is gravita- 
tional clustering, which makes it quite difficult to ex- 
tract any primordial non-Gaussianity from the galaxy 
distribution Durrer et al.(2000), . The CMB on the 
other hand is to first order free of this compli- 
cation and is therefore ideally suited to study the 
early universe. Nevertheless, secondary effects like the 
Sunyaev-Zel'dovich (SZ) effect Aghanim & Forni(1999), 
Cooray(2001), , the Ostriker-Vishniac effect Castro(2002), 
, lensing Cooray & Hu(2001), Bernardeau et al. (2002), 
and others add their own contributions to the total non- 
Gaussianity. To these cosmological effects, we have to add 
foregrounds as well as systematic effects and instrumental 
noise. In order to disentangle all those sources from one 
another, it is essential to assemble a "tool-box" of well- 
understood, fast and robust methods for probing different 
aspects of the non-Gaussian signatures. 

In this paper, we compare the behaviour of two dif- 
ferent classes of estimators, namely the higher order mo- 
ments in wavelet space, and the bi- and trispectrum. Both 
of these do not act directly on pixels, but are applied in 
a dual space. The purpose of this study is not to present 
novel ways of detecting non-Gaussian signatures, but to 
discuss practical aspects of these different approaches, and 
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to compare their behaviour when applied to a selection of 
synthetic "benchmark" maps which are designed to mimic 
typical non-Gaussian contributions to the CMB. 

After the seminal work of Pando et al.(1998) which 
consisted of applying wavelet techniques to look for 
non-Gaussian signatures in the COBE-DMR data, sev- 
eral methods based on the wavelet analysis have 
been developed in the context of the CMB detec- 
tion of non-Gaussian signatures Aghanim & Forni(1999), 
Forni & Aghanim(1999), Hobson ct al.(1999), 

Barreiro & Hobson(2001), Cayon et al.(2001), 

Jewell(2001), Martinez- Gonzalez et al.(2002), 

Starck et al.(2003), . These methods have proven partic- 
ularly suitable for statistical studies due to the combina- 
tion of two characteristics. First, wavelets are quite lo- 
calised both in space and in frequency, which allows for 
the features of interest in an image to be present at differ- 
ent scales. Second, the linear transformation properties of 
Gaussian variables preserve Gaussianity. Conversely, any 
non-Gaussian signal will exhibit a non-Gaussian distribu- 
tion of its wavelet coefficients. The two properties com- 
bined together allow us to associate the statistical signa- 
tures with the spatial features that have caused them. 

The wavelet based estimators have been tested and 
qualified in terms of their sensitivity to the non-Gaussian 
signatures on a variety of simulated data sets such as 
cosmic strings, the SZ effect from galaxy clusters, and 
anisotropics from inhomogeneous reionisation (see the 
above mentioned references). Also the galactic contam- 
ination from dust emission can induce foreground non- 
Gaussian signatures. Wavelets have been used in such a 
context by Jewell(2001) to quantify the predicted level of 
non-Gaussian features in the IRAS maps. Additionally, 
the wavelet based non-Gaussian analysis has also been 
appHed to the COBE-DMR maps. 

Different decomposition schemes and wavelet bases can 
be used for the non-Gaussian studies. In the following, 
we use a bi-orthogonal wavelet basis and a dyadic de- 
composition scheme. This choice has been first motivated 
in Forni & Aghanim(1999) and Aghanim & Forni(1999) 
by the fact that it is the best scheme in the context of 
statistical analysis of CMB signals. It is indeed optimal 
for statistics since it gives, at each scale, the maximum 
number of significant coefficients. In addition, it natu- 
rally allows us to benefit from spatial correlations in the 
signal at each scale. The better performances of the bi- 
orthogonal dyadic wavelet decomposition were confirmed 
by Barreiro & IIobson(2001); and also more recently by 
Starck et al.(2003), in a comparison with the "a trous" 
and ridgelet decompositions. 

We now focus on the Fourier based analysis. The 
interest in using higher order correlation functions to 
study non-Gaussianity dates back many years Luo(1994a), 
Jaffe(1994), Heavens(1998), . The first detection of a non- 
Gaussian contribution in real CMB data was achieved us- 
ing the bispectrum Ferreira et al.(1998), , but later shown 



to be due to a weak systematic effect Banday et al.(2000). 

Numerous works deal with the construction of 
ideal higher order statistical estimators both in spher- 
ical and flat space as well as with their applica- 
tion to real data (CMB and large-scale structure sur- 
vey data). In the CMB context, research focused 
on the two first higher order correlation functions 
in Fourier space, namely the bispectrum (see recent 
work by Santos et al.(2002) and citations therein), and 
its 4-point equivalent, the trispectrum (see IIu(2001), 
Kunz et al.(2001), Cooray & Kesden(2002)). An impor- 
tant advantage of n-point functions over other non- 
Gaussian statistics such as wavelets or Minkowski func- 
tionals is that they are easier to predict for inflation but 
also for secondary sources. This makes them a very pow- 
erful test, explaining their interest for the scientific com- 
munity. 

Indeed, the Boltzmann equation, which describes the 
time evolution of the perturbations in the cosmic fluid, 
does not mix different Fourier modes. The structure of lin- 
ear perturbations is therefore generally simple in Fourier 
space. Furthermore, the projection of the CMB photons 
onto the sky sphere is simple as well, with the square of 
the spherical Bessel functions coupling the Fourier mode 
k and the angular scale, or more directly the multipole £, 
on the sphere. The success of studying the angular power 
spectrum Ce lies in that much of the initial conditions is 
preserved throughout the evolution. Higher order correla- 
tion functions can be treated in much the same way and 
allow for the use of already existing techniques. It is there- 
fore feasible to predict theoretically the expected contri- 
butions to non-Gaussian signatures from different sources. 
This allows us on the one hand to place limits on sources 
of non- Gaussian features if none is detected, and on the 
other hand might enable us to solve the inverse problem, 
namely identifying the origin of the non-Gaussianity, if 
any is found. 

In this work, we concentrate on the flat space ap- 
proach. We use a normalised bispectrum probing all the 
possible triangle configurations which exist in a homoge- 
neous and isotropic universe, and a normalised diagonal as 
well as nearly-diagonal trispectrum which depends on the 
side and the diagonal of the (nearly) trapezoidal configu- 
ration. These are supposed to give us information on the 
scale dependence of any non-Gaussian signal, but contrary 
to the wavelets, no spatial position information is retained. 

We will start in the next section by presenting the data 
sets used for the comparison. In section 3, we examine the 
estimators of the non-Gaussian signatures. The statistical 
analysis and the results are presented in section 4. We 
then end with our main conclusions. 

2. Data sets 

We test the statistical methods to detect non-Gaussian 
signals on three different data sets chosen so that they are 
representative of many astrophysical situations. It is be- 
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yond the scope of the present study to characterise the as- 
trophysical processes represented by the maps. The three 
sets have rather to be considered as test cases used in or- 
der to insure that the non-Gaussian detection is neither 
specific to one pecuhar type of non-Gaussian signatures, 
nor due to one particular detection method. Additionally, 
the combination of the different detection techniques on 
the different sets of data allows us to compare the relative 
sensitivities of the methods on an objective basis. 

In practise, many of the astrophysical situations can be 
roughly classified into the following three different types 
of signals: i) spherically shaped structures standing for 
compact astrophysical sources or beam-convolved point 
sources (e.g. galaxies, galaxy clusters) that we refer to as 
the point sources, ii) strongly anisotropic structures that 
represent filaments and elongated structures (e.g. inter- 
stellar clouds) that we refer to as the filaments, and iii) 
non-linear couplings leading to non-Gaussian signatures, 
referred to as the (or ^-squared) map, as we use the 
superposition of a Gaussian random field and its square. 

In order to characterise the non-Gaussian signatures 
exhibited by all these signals, the maps are compared 
to a set of 99 Gaussian realisations having the same 
power spectrum, referred to as the Gaussian reference set. 
Therefore, the only difference between Gaussian and non- 
Gaussian sets is due to the statistical nature of the two 
processes. All the maps (Gaussian and non-Gaussian) are 
composed of 500 x 500 pixels of 1.5 arc minutes aside. 

2.1. Filaments 

The set of maps for the filamentary structures is a selec- 
tion of 28 IRAS maps which represent the galactic dust 
emission at high latitudes. These maps were chosen on the 
basis of their visual aspect so that they exhibit large and 
bright filaments (Fig. 1, upper middle panel). Therefore, 
they all have a low amount of point sources and they ex- 
hibit elongated structures. In order to focus on the sta- 
tistical properties of the maps, the selection has also been 
made on the basis of their power spectrum. We choose the 
maps showing power spectra with similar shapes regard- 
less of their amplitudes. We then modify the initial IRAS 
maps so that their power spectra are rescaled. Figure 2, 
lower panel, shows the power spectrum of the filaments. 
Thanks to this rescaling, we can consider the 28 modified 
IRAS maps as statistical realisations of the same highly 
non-Gaussian process. This data set has a skewness (third 
moment) of 0.89 ± 0.34 (1 a), the excess kurtosis (fourth 
moment) is 2.63±2.42. The Gaussian reference maps (one 
example is given Fig. 1, lower middle panel) have a skew- 
ness of — 0.01±0.17 and an excess kurtosis of — 0.11±0.19. 
The large fluctuations of the skewness and excess kurtosis 
(also present in the Gaussian maps) arise due to the high 
power on large scales associated with the large structures 
in the IRAS maps. Hence, even the maps produced using 
a Gaussian random process could be easily mistaken as 
being non-Gaussian, proving the necessity of comparing 




Fig. 1. Representative maps of the non-Gaussian signals and 
one of their associated Gaussian realisations with the same 
power spectrum. Upper and lower left panels represent respec- 
tively a point source map and one Gaussian realisation. The 
upper and lower middle panels are for the filaments and an as- 
sociated Gaussian realisation. The upper and lower right panels 
represent the map and a Gaussian field with the same power 
spectrum. 

a given test map with a set of Gaussian reference maps 
with the same power spectrum, rather than relying on a 
theoretical prediction of the Gaussian bounds. For a dis- 
cussion of the Gaussian reference maps, sec section 2.4. 



2.2. Point source maps 

The second data set is referred to as the point sources. It is 
constituted of 50 simulated maps which consist of a distri- 
bution of Gaussian shaped sources having different sizes 
and amplitudes. They are randomly distributed on the 
map and have positive and negative signs (see Fig. 1, left 
upper panel). This signal is aimed at reproducing the typ- 
ical situation of a crowded field of 12.5 square degrees con- 
taining ~ lO'' sources. The field therefore exhibits obvious 
confusion effects. All the simulated maps have the same 
bell-shaped power spectrum (Fig. 2, upper panel) and ex- 
hibit a non-Gaussian signal. This non-Gaussian data set 
has a skewness of —0.160 ± 0.078 and an excess kurto- 
sis of 2.19 ± 0.32. A Gaussian reference set of maps with 
the same power spectrum is simulated (Sect. 2.4, Fig. 1, 
left lower panel). The corresponding Gaussian maps have 
a skewness of —0.003 ± 0.022 and an excess kurtosis of 
—0.005 ± 0.036. The point source maps are hence clearly 
marked as being non-Gaussian by their excess of kurtosis. 



2.3. x-^Quared maps 

The two previously described types of non-Gaussian 
signals are ideal to model astrophysical induced non- 
linearities due to foregrounds and/or secondary effects. 
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Fig. 2. The power spectra of the studied signals in arbitrary 
units. The upper panel shows the power spectrum of both the 
point sources and the maps. The lower panel represents the 
power spectrum of the filaments. 



However in a wider context, we can expect the exis- 
tence of primordial sources of non-Gaussian signatures. 
A way of modelling some of these signatures is by in- 
troducing a non-linear coupling in an originally Gaussian 
distributed perturbation field Coles & Barrow(1987), 
Komatsu & Spergel(2000), Bartolo et al.(2002), . The 
simplest weak non-linear coupling is given by x(^) — 
XlW + /wL(xi(x)- < xiW >), where xl denotes the 
original linear Gaussian field. The non-Gaussian signature 
introduced by the coupling is quantified by the so-called 
non- linear coupling constant /at/, (note that < >= 
0). 

To generate these maps we use a method similar to 
the one followed to create a Gaussian distributed field (see 
section 2.4). The non-Gaussian character is introduced 
by adding to the white noise field its squared minus its 
squared average weighted by /tvl. When multiplying the 
non-linear field in Fourier space with the desired power 
spectrum, a normalisation by 1 -I- 2/^^ is needed in or- 
der to recover the power spectrum from the final map. In 
our approach, Jnl quantifies the amount of the type 
contribution in the Gaussian map. 

We compute two sets of 50 maps for /at^ = 0.01 
and 0.1 in order to test the sensitivity of the statisti- 
cal methods to the amount of non-Gaussian signal intro- 
duced in the Gaussian field. For illustration purposes, we 
show in Fig. 1 right upper and lower panels a map for 
Inl = 0.01 and one associated Gaussian reahsation. The 
power spectrum of these maps was chosen arbitrarily to be 
the same as for the point source maps. As mentioned be- 
fore, these maps are not meant to reproduce accurately the 
physical processes - to simulate a "real" CMB map with 
a primordial x^ contribution we would need to integrate 
over the radiation transfer function. Nevertheless, to relate 
the order of magnitude of our numbers to the ones found 
in the literature, we note that they are rescaled by ap- 



proximately a factor of <i> AT/T ~ 10~^ so that a value 
of /jvL = 0.01 here corresponds to about Jnl ~ 1000 
elsewhere (in reality a bit less, since e.g. on large scales 
AT/T « $/3). The x-squared maps with Jnl = 0.01 
have a skewness of 0.009 ± 0.019 and an excess kurtosis of 
0.002 ±0.031, both of which arc consistent with zero. The 
maps with f^^^ = 0.1 on the other hand have a skewness 
of 0.140 ± 0.024 and an excess kurtosis of 0.037 ± 0.036. 
The Gaussian reference set is the same as for the point 
source maps, see section 2.2. 

2.4. Gaussian realisations 

In the process of investigating the statistical character 
of non-Gaussian signals, as we always deal with a fi- 
nite ensemble of maps, there can be substantial fluctu- 
ations in the results (see for example the IRAS case, sec- 
tion 2.1). Furthermore, although not important in our 
case but rather when analysing experimental data, of- 
ten a "Gaussian" map is not created directly from a pure 
Gaussian random process, but has super-imposed known 
systematic effects (which would result in spurious non- 
Gaussian features). In general, the output of an estimator 
applied to a test-map is therefore not compared to a the- 
oretical result but to the corresponding output obtained 
from a set of Gaussian reference maps (with any known 
systematic effects one may need to add). Any detection of 
non-Gaussian signatures can then be reliably quantified 
by comparing these two sets of results, and we discuss 
this procedure in section 4.1 in more detail. 

We generate two sets of Gaussian distributed fields 
with the power spectrum of the corresponding non- 
Gaussian set. The first one is the aforementioned reference 
set. The second one (called the Gaussian counterpart set) 
is strictly speaking redundant, but is used to illustrate the 
fluctuations arising from finite ensembles, and the pres- 
ence of low-probability results which appear even when 
comparing Gaussian maps with each other if many differ- 
ent tests are applied. 

The simplest standard method Peacock(1999), is to 
create a spatial array of white noise using a random num- 
ber generator with zero mean and unit variance. We can 
then give this Gaussian field the desired power spectrum 
by multiplying it in Fourier space with the square root 
of that power spectrum, the one of the non-Gaussian set 
in our case. Using this procedure, we compute one set of 
149 Gaussian maps (reference set and counterparts) with 
the same power spectrum as the 50 point-source maps, 
one set of 127 Gaussian maps (reference set and coun- 
terparts) with the same power spectrum as the 28 IRAS 
modified maps. We arbitrarily choose the power spectrum 
of the point source maps for the x^ maps and hence use 
the same set of Gaussian realisations. We show in Fig. 
1 (lower panels) one representative Gaussian realisation 
for each data set with the same spectrum as the non- 
Gaussian map. The skewness and the excess kurtosis of 
all the Gaussian reference sets and counterparts are ex- 
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pected to be zero. However, in practice they have a finite 
scatter. The lower limit for the standard deviation of the 
skewness is given for a normal distribution with unit vari- 
ance by yJlh/N'^, or about 0.0077 in the case of N'^ = 500^ 
independent values per map. The expected standard de- 
viation for the excess kurtosis, in the same idealised case, 
is y/OQ/N^ or 0.020. One should note however that the 
scatter of higher order moments is much larger in realistic 
cases and is strongly susceptible to the details of the ran- 
dom process (see sections 2.1, 2.2 and 2.3 for the actual 
numbers for the Gaussian reference sets) . 

3. Estimators of the non-Gaussian signatures 

We examine two large families of estimators for non- 
Gaussian signatures. Both are based on dual space anal- 
ysis. On the one hand, we focus on methods which are 
completely defined in Fourier space (namely the bi- and 
trispectrum). On the other hand, we study estimators in 
the wavelet space (namely third and fourth moments of 
the wavelet coefficients) in which signals are located both 
in pixel and in frequency space. It is worth noting that 
the skewness (third moment) and excess kurtosis (fourth 
moment) of the wavelet coefficients are quite similar to 
respectively the bi- and trispectrum, in Fourier space. 
Additionally, we compute the cumulative probability func- 
tion (CPF) of the pixel distribution. This latter is also one 
of the Minkowski functionals, often called A, the surface 
of the excursion set at a given amplitude. As the CPF is a 
less powerful non-Gaussian estimator, in the sense that it 
does not give information on the scales nor on the physi- 
cal location of any detection, we intend to use it only as 
a baseline for the detection of non-Gaussian signatures on 
a given map. 

As discussed above, we analyse a certain number TVng 
of non-Gaussian as well as Gaussian maps with the same 
power spectrum and compare their properties to a refer- 
ence set of 99 Gaussian realisations, also with the same 
power spectrum. This section describes the statistical es- 
timators we use for our analysis. For a discussion of the 
tools used to interpret the results, see section 4.1. 

3.1. Pixel distribution function in real space 

As we have seen, in the Gaussian case each pixel is an 
independent realisation of a Gaussian random process. 
Therefore, one of the easiest and most direct tests of the 
Gaussian hypothesis just estimates the cumulative prob- 
ability function (CPF) from the ensemble of pixels iVtot- 
To this end, we set 

P(< x) = 7V[pixels with values < x\/Ntot- (1) 

Since we have several independent maps for both the 
Gaussian as well as the non-Gaussian processes, we can 
improve our estimate of the CPF, and hence the sensitiv- 
ity of the test, by just combining the pixels from all maps. 
The two CPFs (from the Gaussian and non-Gaussian pro- 



cesses) can then be compared with statistical methods de- 
tailed in section 4.1. 

Power spectra which vary extremely rapidly (like the 
ones used here) can strongly amplify the random fluctua- 
tions. In this case, the KS test will lead to spurious detec- 
tions of non-Gaussianity even for Gaussian maps (see e.g. 
Aghanim & Forni(1999)). It is therefore best to "whiten" 
the maps first, by deconvolving them with their estimated 
power spectrum. 

Due to its simplicity and its inherent locality (the CPF 
does not test the relative position or configuration of the 
fluctuations, as the location of the pixels in the map is 
not taken into account), and as stated previously, we use 
the CPF method as a baseline for the comparison between 
the non-Gaussian methods. It should be stressed clearly 
that this test can only detect non-Gaussian signatures but 
does not give any further information about their location 
or scale dependence. It can therefore only very weakly 
distinguish between different non-Gaussian contributions. 

3.2. Wavelet space analysis 

The principle 
behind the wavelet transform Grossman & Morlet(1984), 
Daubcchies(1988), Mallat(1989), is to hierarchically de- 
compose an input signal into a series of successively lower 
resolution reference signals and associated detail signals. 
At each decomposition level, j, the reference signal has a 
resolution reduced by a factor of 2^ with respect to the 
original signal. Together with its respective detail signal, 
each scale contains the information needed to reconstruct 
the reference signal at the next higher resolution level. 

The wavelet analysis can be considered as a series of 
bandpass filters. It can thus be viewed as the decomposi- 
tion of the signal in a set of independent, spatially oriented 
frequency channels. Using the orthogonality properties, a 
function in this decomposition can be completely charac- 
terised by the wavelet basis and the wavelet coefficients in 
this decomposition. 

The multi-level wavelet transform (analysis stage) de- 
composes the signal into sets of different frequency bands 
by iterative application of a pair of Quadrature Mirror 
Filters (QMF). A scaling function and a wavelet function 
are associated with this analysis filter bank. The continu- 
ous scaling function (f)A{x) satisfies the following two-scale 
equation: 

0A(a^) = \/2X]^o(ri)0A(2a;-n), (2) 

n 

where /iq is the low-pass QMF. The continuous wavelet 
iPa{x) is defined in terms of the scaling function and the 
high-pass QMF hi through: 

ijA{x) = V2j2hi{n)(pAi.2x-n). (3) 

n 

The same relations apply for the inverse transform 
(synthesis stage) but, generally, different scaling and 
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wavelet functions (0s (x) and -ipg^x)) are associated with 
this stage: 

0s(a;) =:^/2^5oW0s(2x-7^), (4) 

n 

V'sW = \/2^giH0s(2a;-n). (5) 

n 

Equations 2 and 4 converge to compactly supported 
basis functions when 

5]/^o(n) = ^.goH = ^/2. (6) 

n n 

The system is said to be bi-orthogonal if the following 
conditions are satisfied: 

I dpj^{x)cps{x - k)dx = 5{k) (7) 

JR 

[ Mx)Mx - k)dx ^ (8) 

JR 

/ (l)s{x)'ipAix ~ k)dx ^ (9) 
Jr 

Cohen et al.(1990) and Vetterli k Herley(1992) give a 
complete treatment of the relationship between the filter 
coefficients and the scaling functions. 

The wavelet functions are quite localised in space, 
and simultaneously they are also quite localised in fre- 
quency. Therefore, this approach is an elegant and pow- 
erful tool for image analysis, because the features of in- 
terest in an image are present at different character- 
istic scales. Consequently, different wavelet transforms 
have been studied such as the " a trous" algorithm 
Starck et al.(1998), that decomposes a. N x N image / 
as a superposition of the form 

J 

I{x,y) ^ cj{x,y) +^Wj{x,y), (10) 

where cj is a coarse or smooth version of the original image 
I and Wj represents the details of / at scale , and the 
dyadic wavelet transform Mallat(1998), that decomposes 
the signal s in a series of the form : 

,7 

s(0 =5Z^^.^(<^A)j,;(fc)+^^(V'Akz(/eH-,fe (11) 

fc k j=i 

where J is the number of decomposition levels, Wj^k the 
wavelet (or detail) coefficients at position k and scale j 
(the indexing is such that j ~ 1 corresponds to the finest 
scale, i.e. highest frequencies), and cj is a coarse or smooth 
version of the original signal s. 

In the present study and follow- 
ing Forni & Aghanim(1999), we choose among all the pos- 
sible bases and decomposition schemes the bi-orthogonal 
wavelet basis and the dyadic decomposition. We have cho- 
sen to perform a six level ( J = 6) dyadic decomposition 
of our data. In practice, a dyadic decomposition refers to 



a transform in which only the reference sub-band is de- 
composed at each level. In this case, the analysis stage is 
applied in both directions of the image at each decompo- 
sition level. The total number of sub-bands after J levels 
of decomposition is then 3J + 1. This decomposition gives 
as a result the details for the studied signal in both direc- 
tions. At each level, we therefore end up with 3 sub-bands 
that we refer to as the horizontal, vertical and diagonal 
details, plus the smoothed signal. This kind of decompo- 
sition furthermore allows us to benefit from correlations 
between the two directions at each level, and also from the 
maximum number of coefficients, which is crucial for sta- 
tistical analysis. The linear transformation properties of 
Gaussian variables preserve the statistical character mak- 
ing the wavelet coefficients of a Gaussian process being 
Gaussian distributed. Conversely, any non-Gaussian sig- 
nal will exhibit a non-Gaussian distribution of its wavelet 
coefficients. In addition, the multi-scale wavelet analysis 
allows us to associate the statistical signatures with the 
spatial features that have caused them. However, we do 
not investigate the latter property in the present study. 

3.3. Fourier space analysis 

We now describe in some detail the Fourier space methods 
chosen for this comparison. In Fourier space, scalar func- 
tions on the sky sphere, such as the CMB temperature 
fluctuations, are usually expressed as the coefficients aim 
of an expansion in spherical harmonics. For small patches, 
we can instead approximate the sphere by a flat surface. 
This makes it possible to use a fast Fourier transform to 
compute the coefficients. In this case, the transformations 
between pixel space (the original map r(x)) and Fourier 
space are given by 

a(k) = /'d22;r(x)e2'^''*'^ ^ ^r(x)e2"*'^/^'(Ax)2(12) 

r(x) = J d\a{k)e~^'''^'' 

-^l^Y.^{\.)e-'^^^-/^\Akr. (13) 

k 

N'^ is the number of pixels per map. The flat space pixel 
size is determined from the angular pixel size 6 (given in 
arc minutes) via Ax = 27r6'/(360 x 60). The pixel size in 
Fourier space is then Ak ~ 1/ Ax and the angular scale in 
Fourier space is ^ = 27rfc. 

The power spectrum of the temperature distribution 
is the two point function, 

C{k) a(k)a(-k) >w J dipa{k,ip) a{k,~ip) (14) 

The last correspondence is due to the statistical isotropy 
of the temperature field on the sky sphere, which allows 
us to replace the ensemble average with an average over 
directions for a given mode. 

For a Gaussian random field, all information is con- 
tained in the power spectrum. All higher order rt-point 
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functions can be derived from it. All functions with n odd 
vanish, and the even n ones are found via Wick expansion. 
As an example, the four point function is thus found to 
be {a{ki)a{k2)a{k'i)a{k4)) = {a{ki)a{k2)) {a{k'i)a{k4)) + 
{a(ki)a{k^)) {a{k2)a{k4)) + {a{ki)a{k4)){a{k2)a{ky,)) , and 
so on. 

The direct computation of higher order functions in > 
2) for large maps can be very slow. It generically scales as 
^puc^- spherical faster approach has been 

developed over the last years Spergel & Goldberg(1999), 
Hu(2001), De Troia et al.(2003), , by creating maps which 
contain only one scale ^ each. We can use the same ap- 
proach also in flat space, defining the scale maps as 

= ^Ew^^(|k|)«(k)e-'^'^''^/^^Afc)^ (15) 

k 

The window functions Wi serve to select the scale, gen- 
erally we will use Wi{k) = 1 if fc is in the band i and 
otherwise. Other choices are clearly possible. 

Although the formulation of the scale map method is 
more elegant in spherical space, the flat space approach 
has the advantage of using only sums, multiplications and 
fast Fourier transforms of the map, rendering it very quick 
and efhcient. 

With this choice of scale maps, we can compute the 
power spectrum via 

Er£(x)2(Aa;)2 ~ f dkdipkWi{k)a{k,ip)a{k,-ip) (16) 

which contains an additional weight k compared to equa- 
tion (14). In the spherical case, we find instead a Wigner 
3J symbol which has to be taken care of by the nor- 
malisation. This corresponds in the flat space case to 
an infinitely narrow bandwidth (i.e. window functions 
Wg oc S{2Trk — i)). While a finite bandwidth will introduce 
a bias in the result, for sufficiently narrow bandwidths k is 
approximately constant within the band and just provides 
an offset which can also be taken care of via an appropriate 
normalisation; in practice we just divide by the number of 
pixels per band. Tests show that for a bandwidth of two 
pixels we do not see any differences from eq. (14) even 
for rapidly varying power spectra. Similar normalisations 
will have to be performed for the higher order correlation 
functions. 

The bispectrum 

The bispectrum is the three point function in Fourier 
space. In the flat space approximation, it is given by a 
triangle configuration, 

Bg.e.e, =< a(ki)a(k2)a(k3) > 6{ki + ks + kg) (17) 

where the i5 function arises due to statistical isotropy (and 
all other configurations should vanish). 

The above mentioned method of scale maps can be 
directly extended to higher moments. For the bispectrum. 



we compute 

ET,,(x)r,,(x)r,3(x)(Aa.)2 = 

X 

^ E a(ki)a(k2)a(k3) X (18) 

kik2k3 

Wedki)We,{k2)We,{k3){AkfS{k,+k2 + k3). 

Here, we normalised the S symbol so that S{k = 0) = 
l/(Afc)2. 

Again we use the approximation of narrow bands to 
derive the normalisation. We assume that the bispectrum 
does not vary strongly within a given band. Then the nor- 
malisation which needs to be divided out is given by 

N{h,£2, 4) = Yl i\q\)Wi, (|k + q|), (19) 

k,q 

which we sum up using our choice of scale window func- 
tions. Since this function is geometrical in nature and in- 
dependent of the actual map in question, it needs to be 
computed only once. This is again just the number of con- 
figurations which contribute to a given (£1,^2,^3) triplet. 

The variance of this estimator for a Gaussian random 
field is then found to be cr| = 6Ci,Ci^CejN{£i, £2, £3) 
(for our Wi which verify W^{k) = W,(k)). 

We will in general divide out the dependence of the 
bispectrum on the power spectrum, and we will do the 
same for the trispectrum. This Ci normalised estimator, 
Bi^i^g^/lCe^Ce^Ci^}^^"^ has comparable properties to the 
unnormalised one, but is more robust with respect to 
wrong estimations of the power spectrum. Komatsu et al. 
(2002) even found it to have a slightly lower variance. 

The trispectrum 

The trispectrum, the four point function, is slightly more 
complicated. Statistical isotropy implies again a (5-function 
over the four momenta, so that we can visualise it as a 
quadrilateral with sides £1, £2, £3 and £4, and we can ad- 
ditionally specify (within limits) the length of one of the 
diagonals, a. This additional degree of freedom has to be 
taken into account when constructing an estimator. To 
this end, we fix a diagonal a with length |a| = a. We then 
decompose the full trispectrum into two triangles, aver- 
aged over all directions of a: 

TiMi. = E ri-'^VFa(|a|)(Afc)2+permut.(20) 

a 

i - £ ■ 

The sub-estimators Ta" ^ need to be constructed so 
that one side is given by a and that the other sides have the 
lengths indicated, e.g. £i and £j. This defines the triangle 
completely. A possibility for r^^'^^ is to use the Fourier 
transform of ^ • Tg^ , 

=Y,T,,{^)T,^{^)e^^^'^^/^\Axf . (21) 
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Then, the integral in the definition of the Ti (see eq. (15)) 
over the exponentials leads to S{ki + k2 — a)(5(k3 + k4 + a) 
and ensures the correct geometric structure. 

In this study, we are not going to compute all possible 
trispectra. Instead, we are going to concentrate on two 
simpler cases. The first one is the case £i ~ ^2 — ^3 = ^4, 

rp{0) rpa (00\ 

This is quite a natural choice, and due to its symmetry 
well suited to investigate elongated structures like fila- 
ments. Furthermore, one needs only one scale at a given 
time and does therefore not need to keep all scale maps in 
memory, greatly decreasing the amount of memory needed 
by the algorithm. But this "diagonal" estimator docs not 
vanish in the case of a Gaussian random field. It is in 
general proportional to the square of the power spectrum. 
Specifically, 

{t^°I)^ (X C| {N{£fSa,o + 2N{£, a)) . (23) 

Where the normalisation function for our choice of geom- 
etry is given by 

N{£, a) = ^ W^a(|a|)I^,(|k|)VK,(|a + k|), (24) 

k,a 

and again needs to be calculated only once. For a = 
it is N{£,a = 0) = N{£) = Ek^^^d^l), the number of 
non-zero points in the scale map £ (for our choice of scale 
window functions W). This could hide a real non-Gaussian 
contribution behind the Gaussian "noise" . 

It is also noteworthy that in this simple case the di- 
agonal is bounded by < a < 2^. As the length of the 
other diagonal b is related to a via b = ^ {M^ — a^)-, the 
cases with a > \f2l are redundant. But since the binning 
in a and b is different, there is no easy way of actually re- 
moving these cases. The canonical scale decomposition in 
spherical space allows the orthonormalisation procedure 
described below to remove the superfluous values in that 
case. 

We also use the near-diagonal case 

which vanishes (for a 7^ 0) for a Gaussian random field. 

As discussed in Kunz et al. (2001), the trispcctrum 
is not unbiased. This does not concern us overly much, 
since we only want to detect non-Gaussianity in a first 
step, not estimate its precise amplitude. Nonetheless, as 
discussed in that paper, the estimator given here has 
many superfluous degrees of freedom which can influ- 
ence the degree of detection, depending on where pre- 
cisely a signal is located in {£, a) space. We can decrease 
the number of estimators (and so the superfluous de- 
grees of freedom) as follows: We define for each £ the 
matrix Mat =< T^T^ >g, where <>g is the expecta- 
tion value over a Gaussian ensemble, and diagonalise it, 
Mab = J2i ^^''^a''^b- Smce many of the eigenvalues A* vanish, 
we define a smaller ensemble of estimators by the linear 
combinations T; = Ea^/^)^" if 7^ 0- 
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Fig. 3. Two distributions of near diagonal trispectrum values 
for the filaments: The upper panel is for {I, a) — (1957,0), a 
highly non-Gaussian case (dashed line) and its corresponding 
values issued from the Gaussian set of maps (solid line). The 
lower panel is for {£, a) — (110, 400), a nearly Gaussian case. In 
this case, the two distributions (same linestyle as upper panel) 
are very close to each other. 



4. Results 

4.1. Statistical analysis 

In the last section, we have described in great detail the 
estimators which we apply to our synthetic data. Before 
we can discuss the results, it is necessary to explain in 
which way we quantify statistically the level of any non- 
Gaussian detection. To this end, let us study each "test" 
separately. A "test" is each hi- or trispectrum configura- 
tion (as defined by the triplet or quadruplet of £ values), 
and the skewncss or excess kurtosis for the coefficients of 
each combination of scale and orientation in the wavelet 
case. Hence, each of our four estimators is made up of 
many of these tests. We apply each test independently to 
the non-Gaussian maps, their Gaussian counterparts as 
well as to the Gaussian reference set. We therefore end up 
with three ensembles of values for each test. 

We illustrate this by plotting in Fig. 3 (upper panel) 
two distributions of near-diagonal trispectrum values 
with {£, a) = (1957, 0) for both the filament maps and the 
corresponding Gaussian reference maps. The lower panel 
shows the two distributions for {£, a) = (110, 400). Clearly 
in the first case the two distributions are very different, 
corresponding to a detection of non-Gaussianity. It is this 
difference that needs to be quantified. 

A standard method to statistically compare two dis- 
tributions is to perform a Kolmogorov-Smirnov (KS) test. 
The KS test returns the distance d which is defined as 
the maximum value of the absolute difference between 
the two cumulative distribution functions. One can de- 
rive the probability Pks that the two data sets are drawn 
from the same distribution (i.e. the two distributions 
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are statistically the same) using the following equation 
Press et al.(1992), : 

PifS = Qks([v/^ +0.12 + 0.11/ v/iVe]d), (26) 
with 



oo 

?Ks(a;) = 2^(-l)^-icxp(-2jV 



'■) 



(27) 



where Nc 



N1N2 



is the effective number of data points, 



N1+N2 

given A'^i and N2 the number elements of the two distribu- 
tions. In our study, we compare the distribution of values 
for a non-Gaussian data set to a distribution obtained 
from the Gaussian reference set. Therefore, the values of 
the KS probabilities Pks represent a measurement of the 
confidence level, given by 1 — Pks, for the detection of the 
non-Gaussian signatures. For example, if for a certain test 
Pks = 10~^ then this test exhibits non-Gaussian features 
at a confidence level of 99.9%. The probabilities Pq foimd 
by comparing the ensemble of Gaussian counterparts with 
the Gaussian reference ensemble are uniformly distributed 
between and 1. This, after all, is the definition of the 
probability Pks- 

If we now move back from a single test to the estima- 
tors as a whole, we are confronted with a new problem 
which is most noticeable in the case of the bispectrum, 
for example. There, we perform of the order of rit ~ 10^ 
tests. Hence, we expect to find some tests (although very 
few) where even the comparison of two Gaussian ensem- 
bles results in a very low probability of order l/nt- This 
is illustrated e.g. by figure 5 where the uniform distribu- 
tion of the Gaussian probabilities is also clearly visible. 
Only a deviation from this behaviour signals a violation 
of the Gaussian hypothesis. As a corollary, if a single test 
(a single scale) has a probability Pks ^ l/^-t, that test 
alone ensures with near-certainty that the process which 
created the map is not Gaussian. To avoid any confusion 
with these small Gaussian probabilities, we state in gen- 
eral both the results of the KS test of non-Gaussian versus 
Gaussian reference maps as well as the Gaussian counter- 
parts versus Gaussian reference results. 

Given the large ensemble of Pks values for the bl- 
and trispectrum, and knowing their distribution in the 
Gaussian case (either taking the theoretical uniform dis- 
tribution or, which is more appropriate, directly the mea- 
sured values), it is tempting to go one step further and to 
apply the KS test a second time. This meta- statistics then 
returns one global meta-probability for the detection of a 
non-Gaussian signature by a given method, but neglects 
many aspects such as correlations between the tests which 
make up that method. In the case of the bi- and trispec- 
trum, the meta- statistics discards the information on the 
frequency location of the tests and thus weakens the power 
of the Fourier analysis. In addition, the meta- statistics 
might also degrade our ability to detect the non-Gaussian 
signatures due to the known limitations of the KS test. 
Even if, as an example, a few locations have very small 



probabilities Pks l/'^ti the global meta-probability may 
be acceptable if those points are not numerous enough. Of 
course these highly significant values alone would suffice 
to reject the Gaussian hypothesis. 

It is worth noting that the choice of the number of 
test (non-Gaussian) and reference maps for the first KS 
test has an important impact on the meta- statistics. If the 
number of maps in the two sets have not enough different 
prime factors (e.g. 50 and 100), then heavy aliasing can 
occur since the KS test is only sensitive to the difference 
in the values of the cumulative distribution functions. In 
this case, only very few distinct differences (and hence KS 
probabilities) exist, leading to a non-continuous distribu- 
tion of the KS probabilities. As a consequence, the second 
KS test gives a spuriously low global meta-probability in- 
dicating that the signals are statistically quite different, 
even when comparing Gaussian versus Gaussian maps. To 
avoid this, we use 50 test maps and 99 reference maps, 
which largely eliminates the problem. 

For the wavelet based estimators, we could also per- 
form meta- statistics. We would in this case have to cre- 
ate a distribution of KS probabilities for each decompo- 
sition scale by comparing several sets (Gaussian versus 
Gaussian and non-Gaussian versus Gaussian), and then 
apply the KS test a second time. We refrain from apply- 
ing this method here. 

The KS test tends to be more sensitive around the me- 
dian value of the cumulative distribution function and less 
sensitive to differences at the extreme ends of the distribu- 
tion. It is therefore a good test to measure shifts but not 
so good in finding spreads which might affect the tails of 
the distribution. A way of accounting for these differences 
at the tails is to replace the KS distance d by a stabilised 
or weighted statistics, for example the Andersen-Darling 
(AD) test or the Kuiper test. In the former test, no prob- 
ability is computed. This property makes it vital for the 
estimation of the non-Gaussian detection level that we 
actually compare on the one hand Gaussian counterparts 
and the Gaussian reference set, to on the other hand non- 
Gaussian data set and the reference set. We have checked 
that both the AD and the Kuiper tests give results that are 
in quite agreement with the KS test. In the following, we 
concentrate onto the KS test, and we present the results 
for each statistical estimator in terms of the probability 
that two distributions are identical. 



4.2. Results from the cumulative probability function 

As mentioned before, our baseline is given by a sim- 
ple comparison between the estimated CPF of the non- 
Gaussian data sets and the corresponding reference set. 
We use the KS test for this comparison, since it provides 
us directly with an estimate of the probability. We do not 
collect all pixels of all test maps, but use only the pixels 
of one map at a time and find therefore a distribution of 
values. Table 1 gives the averaged results obtained in this 
way, compared to the results from the Gaussian counter- 
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Pks 
Pg 



Filaments Point sources 



X maps 




0.73 




0.72 



4.7- 10" 
0.72 



f recovered from white noise maps by the sl<ewness 



Table 1. Test for non-Gaussianity using the CPF for all 
(whitened) maps. The first line is the average KS probability 
for the real, non-Gaussian maps, the second line the average 
probability for the Gaussian counterpart set. The type maps 
have a non-linear coupling factor Jml ~ 0.01. 



part set. The maps used for this purpose were deconvolved 
with the average power spectrum extracted from the com- 
parison set. 

In all cases, the non-Gaussian signature is clearly de- 



tected (with the fNL ~ 0.01 x type maps being the 
least non-Gaussian by far), while the Gaussian counter- 
part maps are consistent with the Gaussian hypothesis. 
As the point source maps and the maps use the same 
Gaussian reference and counterpart maps, they also share 
the same Gaussian probability Pq- That value might be 
slightly higher than expected due to the deconvolution 
which has to use an estimated power spectrum (the true 
one being unknown in practice). But in the worst case, 
we expect this measure to lead to a weaker detection of 
non-Gaussianity. If wc do not deconvolve the maps, we 
find, Pg = 10~^ for the point source Gaussian reference 
set versus the Gaussian counterpart set, and the filament 
case is much worse. 

Additionally, we can quantify the sensitivity of this test 
to the non-Gaussian signatures present in the x^ maps. To 
this end, we form a likelihood estimator C[fNL], and find 
the best- fit parameter /atl by minimising C. The vari- 
ation of the recovered Jnl for a set of Gaussian maps 
gives an indication of the sensitivity of the approach. A 
priori, we could use directly the KS probability as given 
above for the likelihood function. But, as demonstrated in 
Santos et al.(2002), it is simpler to use a x^ fitting proce- 
dure to a signal which can be analytically predicted and 
in which /nl enters linearly. In this case, we can calculate 
the best-fit value in a very straightforward way. 

The probably simplest signals which can be extracted 
from the CPF are its moments. As the average vanishes, 
and as all even moments show only a leading order non- 
Gaussianity proportional to ^ 1, it is best to use the 
skewness. To first order in Jnl, we find < >= Q/nl'^'^, 
where the standard deviation a is the one of the original 
Gaussian map. The one recovered from a non-Gaussian 
X^ type map differs by a term of order which can be 
neglected. The likelihood is then given by 



(r(x)3 - /jv^6a^ 
15ct6 



(28) 



Since it is again important to use a whitened map we will 
assume that a = 1. We set dx^ /dfNL = and find 



1 




3 

X 10"' 



Fig. 4. The values of fNL recovered from the skewness of 1000 
maps which contain Gaussian white noise. The variance is 



A/a 



8. 10" 



(29) 



Figure 4 shows the recovered Jnl from 1000 maps with 
Gaussian white noise. Their variance is cr = 8. 10~^ and 
we find for the 50 maps with /jvl = 0.01 that /nl ~ 
0.0098 ±0.0007. 

There is a caveat: In order to extract that value, we 
needed to know the exact kind of non-Gaussian signal we 
are dealing with. However, it is valid to use the present 
procedure to set limits on the parameter in question if no 
non-Gaussian signal is detected. 

4.3. Results from the wavelet decomposition 

In order to find the non-Gaussian signatures present in 
the studied signals, wc compute the skewness and excess 
kurtosis of the wavelet coefficients at each decomposition 
scale and in each sub-band, i.e for the three types of de- 
tails (horizontal, vertical and diagonal). We compute these 
quantities for both the Gaussian and non-Gaussian maps 
and we compare the obtained distributions of skewness 
and excess kurtosis for the Gaussian and non-Gaussian 
processes using the KS test. This comparison allows us to 
get, for each scale and orientation, the probability that the 
non-Gaussian process and the Gaussian realisations with 
the same power spectra have the same distribution. 

All the obtained probabilities for the skewness and ex- 
cess kurtosis are given in the tables 2, 3 and 4. In these 
tables and for each type of details, the first set of 2 lines 
represents Pks for the skewness. The first line is for non- 
Gaussian versus reference set comparison, and the second 
line is for the Gaussian counterpart versus reference set 
(G vs G). The second set of 2 lines represents the Pks 
for the excess kurtosis. Again, the first line is for non- 
Gaussian versus reference set comparison, and the second 
for the Gaussian counterpart versus reference set. 
The probabilities of the order of, or lower than, 10^^ indi- 
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eating a non-Gaussian detection at a level > 99.9% are all 
in shaded boxes. The probabilities obtained from the com- 
parison between the Gaussian counterparts and the refer- 
ence Gaussian realisations (G vs G) allow us to estimate 
any possible non-significant detection of non-Gaussian sig- 
nature, or statistical fluctuation in the Gaussian sets. 

4.3.1. Third moment of the wavelet coefficients: 
Skewness 

Confirming the CPF results for which there was a clear de- 
tection of non-Gaussian signatures in the set of maps 
with a coupling factor /nl = O.f , the non-Gaussian signa- 
tures are very easily detected in the wavelet space. We 
therefore focus on the maps with a coupling factor 
/nl = 0.01. We summarise the KS probabilities derived 
from the comparison of the distribution of skewnesses in 
Table 2, first two lines of each type of details. For these 
maps, the departure from the Gaussian hypothesis is 
only detected at the first wavelet decomposition scale (i.e 
3 arc minutes) . Both the vertical and horizontal details ap- 
pear very sensitive to the non-Gaussian features present 
in the maps. The non-Gaussian signatures are also de- 
tected in the diagonal details, at the first scale, but with 
a much lower confidence level. From the KS probabilities 
obtained from the comparison Gaussian counterparts ver- 
sus Gaussian reference maps (second line of the first set 
in each detail) , it is obvious that the detection of the non- 
Gaussian signatures are highly significant. It is however 
worth noting that we detect some fortuitous non-Gaussian 
signatures, at 93 to 96% confidence levels, due to statisti- 
cal fluctuations when we compare the Gaussian maps with 
the reference Gaussian set. 

In the case of the filaments (Table 3, first two lines of 
each type of detail) , the KS test for the skewness distribu- 
tion of the non-Gaussian against the Gaussian reference 
set shows that the non-Gaussian features are detected at 
all scales (with less confidence at scale 6). The detection 
takes place for the three types of details. We note that 
the probability for the non-Gaussian signal to be com- 
patible with a Gaussian process increases for the highest 
decomposition levels, i.e. the largest angular scales, as it is 
expected from our studied signal in which the filamentary 
structures have not increasingly large sizes. 

As for the point sources (Table 4, first two lines of each 
type of details) , the non-Gaussian signatures do not show 
up very significantly (whatever decomposition scale and 
details). The obtained probabilities remain always larger 
than ~ 0.002 suggesting a non-Gaussian detection at a 
~ 99.8% confidence level at best. 

For the filamentary structures as well as for the point 
sources, the comparison between Gaussian counterparts 
and Gaussian reference (G vs G) set shows that we are 
very much compatible with identity test (processes hav- 
ing the same distribution). The obtained KS probabilities 



show the high level of significance of the non-Gaussian 
detection. 

4.3.2. Fourth moment of the wavelet coefficients: 
Excess kurtosis 

We now turn to the fourth moment in the wavelet space, 

1. e. the excess kurtosis of the wavelet coefficients. The re- 
sults are given by the second sets of two lines in Tables 

2, 3 and 4. Again the first line is for the KS probability 
obtained from the comparison of non-Gaussian maps with 
Gaussian reference set, and the second line stands for the 
comparison of Gaussian counterparts versus reference set 
(G vs G). 

In the case of the x^ maps with /nl = 0.01, the 
fourth moment do not, as expected, exhibit any depar- 
ture from the Gaussian hypothesis. This is demonstrated 
by the large values of the KS probabilities (Table 2, second 
sets of lines) . 

The point sources (Table 4, first line of the second set 
of lines for each detail), exhibit very highly non-Gaussian 
signatures as shown by the small KS probabilities (always 
smaller than ~ 10^^). Furthermore, the non-Gaussian sig- 
natures are present at all decomposition scales and for all 
the details. We compare the obtained probabilities with 
those derived from the comparison of Gaussian counter- 
parts with Gaussian reference set. We note the high level 
of significance of the non-Gaussian detections. 

The same kind of results (detection of non-Gaussian 
features at all scales with all details) are obtained for the 
filaments (Table 3, second sets of lines) for which the prob- 
ability that the non-Gaussian maps are statistically equiv- 
alent to the Gaussian reference maps is always smaller 
than 10^^. We have again computed the KS probabili- 
ties for the excess kurtosis of the wavelet coefficients from 
the comparison of Gaussian maps. This comparison shows 
that the two sets of maps are statistically the same which 
conversely reinforces the level of significance of the non- 
Gaussian detection in the filament maps. 

4.4. Results from the bispectrum 

Similarly to the wavelet coefficients, we compare for each 
triplet (^1, ^2, ^3) the distributions of the (normalised) bis- 
pectrum values, as described in section 3.3, for the non- 
Gaussian realisations and for the Gaussian reference set. 
Hence, we associate to each triplet {^1,^2^3) two ensem- 
bles of bispectrum coefficients - on the one hand 50 or 28 
values for the non-Gaussian maps, and on the the other 
hand 99 values for the Gaussian reference set. Using the 
KS test, we therefore end up for each triplet with a prob- 
ability that the signals are drawn from the same process, 
i.e. that the non-Gaussian maps are compatible with the 
Gaussian reference set. 

For the full bispectrum, we limit ourselves to a band- 
width of four pixels, but probe additionally the diagonal 
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Scale 1 


Scale 2 


Scale 3 


Scale 4 


Scale 5 


Scale 6 


Vertical 


2.0810"^® 


0.226 


0.488 


0.828 


0.158 


0.368 


G vs G 


0.411 


0.450 


0.160 


0.160 


0.791 


0.036 




0.488 


0.426 


0.367 


0.426 


0.426 


0.189 


G vs G 


0.332 


0.791 


0.068 


0.940 


0.411 


0.332 


Horizontal 


1.97 10~" 


0.002 


0.828 


0.555 


0.695 


0.625 


G vs G 


0.876 


0.207 


0.940 


0.876 


0.940 


0.332 




0.315 


0.0450 


0.625 


0.764 


0.426 


0.226 


G vs G 


0.264 


0.596 


0.411 


0.940 


0.596 


0.791 


Diagonal 


0.001 


0.157 


0.828 


0.368 


0.315 


0.315 


G vs G 


0.791 


0.596 


0.940 


0.695 


0.049 


0.791 




0.828 


0.695 


0.930 


0.828 


0.764 


0.963 


G vs G 
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Table 2. For the ^ maps with /nl = 0.01: The KS probability (for two signals to have identical distributions) for each 
of the details at all decomposition scales. For each detail, the first two lines are for the skewness of the wavelet coefficients, 
and the second two lines represent the excess kurtosis. For each pair, the first line stands for the comparison between the 50 
non-Gaussian and 99 reference Gaussian maps, and the second line for the comparison between the 50 Gaussian counterparts 
and 99 reference Gaussian maps (G vs G). All probabilities lower than 10"'' are in shaded boxes. 
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Table 3. For the filaments; The KS probability (for two signals to have identical distributions) for each of the details at all 
decomposition scales. For each detail, the first two lines are for the skewness of the wavelet coefficients, and the second two 
lines represent the excess kurtosis. For each pair, the first line stands for the comparison between the 28 non-Gaussian and 
99 reference Gaussian maps, and the second line for the comparison between the 28 Gaussian counterparts and 99 reference 
Gaussian maps (G vs G). All probabilities lower than 10""^ are in shaded boxes. 



elements with a bandwidth of two pixels. We have tested 
that using the higher resolution for the full bispectrum 
does not lead to a better detection of the non-Gaussian 
signal for the type maps, and does not improve the 
variance of the recovered value of /jvl (see the following 
subsection). 



4.4.1. 



X maps 



For the sake of a clearer visual representation of the re- 
sults, we choose to show the KS probabilities as a function 
of the surface and the smallest angle defined by the triplet, 
rather than as a function of the triplet itself. Note that we 
are speaking of surfaces and angles in Fourier space. 



We first focus on the case of the non-Gaussian ^ maps 
with a coupling factor of /nl = 0.01. We find, for the full 
bispectrum, that the probabilities for the non-Gaussian 
signal to be compatible with its Gaussian reference set 
are as low as a few times 10^^. However as mentioned 
in section 4.1 (see also Fig. 5), the comparison of the 50 
Gaussian counterpart realisations to the reference set of 
99 Gaussian maps gives KS probabilities of the same or- 
der, and it is hard to say whether there is a detection of 
non-Gaussianity or not. If we apply the KS test a second 
time, now to the two resulting distributions of probabili- 
ties, the global meta-probability obtained is 4.1 10"*. The 
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Table 4. For the point sources: The KS probability (for two signals to have identical distributions) for each of the details at 
all decomposition scales. For each detail, the first two lines are for the skewness of the wavelet coefficients, and the second two 
lines represent the excess kurtosis. For each pair, the first line stands for the comparison between the 50 non-Gaussian and 
99 reference Gaussian maps, and the second line for the comparison between the 50 Gaussian counterparts and 99 reference 
Gaussian maps (G vs G). All probabilities lower than 10^"^ are in shaded boxes. 
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Fig. 5. For the maps with non-linear coupling factor /jvi = 
0.01: Distribution of the KS probabilities obtained by i) com- 
paring the full normalised bispectrum estimator of the non- 
Gaussian maps to the Gaussian reference set (dashed line), 
and ii) comparing the same quantity for the Gaussian counter- 
parts to the Gaussian reference set (solid line). In the second 
case, the probabilities are distributed uniformly. The slope is 
due to the log-log representation we have chosen. 



bispectrum therefore does indeed detect this type of non- 
Gaussian signatures. 

As the bispectrum shares many properties with the 
skewness, we would expect it to be rather good at detect- 
ing type non-Gaussianity. Given the way in which we 
construct our test maps, the theoretical signal is 



(30) 



Wc can eliminate any dependence on the power spectrum 
by using the Ci- normalised estimator. The likelihood 



function is then (denoting a triplet (^1,^2,^3) by a) very 
similar to the one for the CPF, eq. (28), 



NL\ 



E 



(31) 



and since the dependence of dy;^ /dfNL = on f^L is hn- 
ear we can again solve for it analytically. By using this 
method, wc recover /nl = 0.010 ±0.002 and find that the 
1000 Gaussian white noise maps have a variance of 0.0023 
(see Fig. 6). As mentioned earlier, using a bandwidth of 
two instead of four does not improve the sensitivity. One 
should not forget, though, that the Fourier space methods 
are better adapted to spherical spaces, where the scale 
decomposition is automatic and does not need to be im- 
posed by placing rings on a rectangular grid. Furthermore, 
the skewness is always just a number, while a realistic 
model of primary non-Gaussian signatures (taking into 
account the radiation transfer function) introduces a non- 
trivial shape and even sign changes into eq. (30), see 
e.g. Santos et al.(2002). In the presence of several different 
kinds of non-Gaussian signatures, the bispectrum should 
be better at discriminating between them in that case, 
and also using a lower bandwidth may become useful. 

4.4.2. Point sources 

In the case of the non-Gaussian maps made of a distribu- 
tion of point sources, the KS test gives again probabilities 
down to 10~^ for the full bispectrum, apart from a sin- 
gle point with much lower probability {Pks — 10~^). We 
concentrate onto the full bispectrum tests and we need 
to compare the obtained KS probabilities with the ones 
derived from the comparison of 50 Gaussian counterparts 
versus the reference set. The figure 7 displays the distri- 
bution of KS probabilities. We note in particular that the 
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Fig. 7. For the point sources: Distribution of the KS probabil- 
ities for the full bispectrum. The dashed line represents the re- 
sults from the comparison of non- Gaussian maps to Gaussian 
Fig. 6. The values of Jnl recovered from the bispectrum of reference set. The solid hue is for the comparison Gaussian 
1000 maps which contain Gaussian white noise. The variance counterparts versus reference set 
is A/jvL = 2. 10"^ 



non-Gaussian signal exhibits an excess of low KS proba- 
bilities as compared with the probabilities from Gaussian 
vs reference set. This difference between the two distri- 
butions is confirmed by the meta- statistics which gives a 
probability of 8.6 10"''. 

Let us now concentrate onto the KS probabilities that 
are lower than 10""^. We display in the surface-angle plane 
the results of the non-Gaussian vs reference set compari- 
son (Fig. 8, diamonds) and Gaussian counterparts vs ref- 
erence set (Fig. 8, asterisks). We notice that for small sur- 
faces (in £ space, i.e. Fourier space) the majority of points 
are associated with the non-Gaussian maps. There are no 
detections associated with the Gaussian vs reference set 
comparison. This is even more obvious when we plot the 
distribution of surfaces for KS probabilities < 10"'^ (Fig. 
8 lower panel, histogram of surfaces). The first bin is in- 
deed dominated by the detections from the non-Gaussian 
versus reference set comparison. The small surfaces in 
Fourier space are associated with large surfaces in real 
space {x cx 1/^). The result wc obtain indicates that we 
are actually detecting the non-Gaussian signatures asso- 
ciated with the higher amplitude spots in the maps (Fig. 
1, upper left panel) that are also spatially more extended 
and lying on a crowded background of smaller sources. 
The histogram of angles (Fig. 8, middle panel) in turn 
does not show any highly significant difference between 
the non-Gaussian versus reference and Gaussian counter- 
part versus reference. 

4.4.3. Filaments 

Our last non-Gaussian data set is constituted of the 28 
modified IRAS maps representing the filamentary struc- 
tures. This process is by far the most non-Gaussian. The 



KS test between the distribution of bispectrum values 
for the 28 non-Gaussian maps and 99 Gaussian reference 
maps goes down to 10~^° for the full bispectrum. Similarly 
to the other processes, we also compare 28 Gaussian reali- 
sations to the reference set in order to estimate the signif- 
icance of the detections. In the case of the filaments, there 
is no doubt that the non-Gaussian signatures are detected. 
This is exhibited (Fig. 9) by the difference in the distri- 
bution of KS probabilities obtained when comparing the 
non-Gaussian maps to the Gaussian reference set (dashed 
line), and the probabilities resulting from the comparison 
of 28 Gaussian realisations versus the Gaussian reference 
set (solid line). Unsurprisingly, the meta- statistics returns 
a probability of confirming the very high significance of 
the non-Gaussian detection. 

In Fig. 10, we display both the Gaussian vs Gaussian 
(gray triangles), and non-Gaussian vs Gaussian (black di- 
amonds) points. We note that most of the probabilities 
in the Gaussian vs Gaussian case lie above the 10"'^ limit 
(horizontal solid line). The non-Gaussian character is very 
strongly localised in the triplets with small surfaces. The 
distribution as a function of angle is much more uniform. 

If we study the KS probability as a function of both 
angle and surface (see figure 11), we find that the non- 
Gaussian signatures are mainly localised in two distinct 
"plumes" . The first contains triplets with angles basically 
no larger than 0.4 radian and surfaces that arc as large 
as 10^, while the second, larger plume, covers all angles 
but only triplets with a slightly smaller surface. The two 
histograms, in Fig. 11, show the number of triplets within 
a range of KS probability as a function of either surface 
or smallest angle. They show that the very lowest proba- 
bilities (solid line) are associated with the second, larger 
plume. They are distributed over a wide range of angles, 
but have all a small surface. The intermediate probabilities 
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Fig. 8. For the point sources: Tiie KS probabilities lower than 
10"'^. In the surface-angle plane (upper panel) defined by the 
triplets, the asterisks are for the comparison Gaussian counter- 
parts versus reference set, and the diamonds for the compari- 
son non- Gaussian maps versus reference set. The middle panel 
shows the distribution of angles for the above selected KS prob- 
abilities. The lower panel shows the distribution of surfaces for 
the same selection. The dashed line is for non-Gaussian versus 
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Fig. 9. For the filaments; Distribution of the KS probabilities 
for the full bispectrum. The dashed line represents the results 
by comparing the non-Gaussian maps to the Gaussian refer- 
ence set. The solid line is for the comparison Gaussian coun- 
terparts versus reference set. 
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Fig. 10. The filament maps: KS probabilities from the compar- 
ison of the full bispectrum as a function of the surface defined 
by the triplet. Black diamonds represent the results obtained 
from the comparison non- Gaussian maps versus Gaussian ref- 
erence set. Gray triangles represent the results obtained from 
the comparison Gaussian counterparts versus reference set. 
The solid horizontal line indicates the limit for Pks = 10~^. 



are concentrated in triplets with small angles but a wider 
range of surfaces. These figures show that the bispec- 
trum allows very detailed analysis of the non-Gaussianity 
present in a data set, if it is detected at a sufficiently high 
level. 

4.5. Results from the trispectrum 

Following the same approach as for the bispectrum and 
as explained in detail in section 4.1, wc compare for each 
pair (£, a) the distribution of the trispectrum values (de- 
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Fig. 11. For the full bispectrum of the filaments: Distribution 
of KS probabilities in the surface-angle plane defined by the 
triplets (upper panel). Diamonds are for the comparison non- 
Gaussian maps versus Gaussian reference set. Black diamonds 
are for the KS probabilities comprised between 10~^ and 10""^, 
gray diamonds represent 10~^° < Pks < 10~^, and light gray 
diamonds stand for Pks < 10^^". The asterisks are for the 
comparison Gaussian counterparts versus reference set. The 
middle panel shows the distribution of angles for the above se- 
lected KS probabilities. The lower panel shows the distribution 
of Riirffl.r.es for the sa.me selectinn, Da.Rherl lines are for in-5 ^ 
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Fig. 12. For the maps with non-linear coupling factor 
Jnl ~ 0.01: Distribution of the KS probabilities obtained by 

i) comparing the trispectrum estimators of the non-Gaussian 
maps to those of the Gaussian reference set (dashed line) , and 

ii) comparing the trispectrum estimators of the Gaussian coun- 
terparts to those of the Gaussian reference set (solid line). 



fined as in section 3.3) for the non-Gaussian set and for 
its Gaussian reference set. We compare the sensitivity of 
the two estimators T'"' and T'^+' (respectively diagonal 
and near-diagonal trispcctra) in terms of non-Gaussian 
detection on the point source maps since the signal is non- 
Gaussian but not as extreme as in the filaments. 

As a significance discriminator, we use the equivalent 
results from the comparison between the set of Gaussian 
counterpart maps and the Gaussian reference set. We 
checked that the normalisation of the trispectrum to the 
power spectrum does not affect the resulting KS probabil- 
ities significantly. 

4.5.1. maps 

For the /nl = 0.01 maps and Gaussian maps, the 
smallest KS probabilities are a few times 10~^, see Fig. 12. 
The meta- statistics confirms that we do not sec any sig- 
nificant deviation from the Gaussian hypothesis: The or- 
thonormaliscd T^^^ estimator leads to a meta-PKS = 0.13 
for a non-linear coupling constant of 0.01. 

It is not surprising that the trispectrum does not de- 
tect any non-Gaussian signature. Just as for the excess 
kurtosis of the pixel distribution in real space (see section 
4.2), the lowest order contribution to the trispectrum is 
proportional to /^^. This makes it naturally much less 
suited to find this kind of non-Gaussian signals. This fact 
can however be used as a counter-check, since if we be- 
lieve to have detected a type signal in the bispectrum, 
we should not find it when applying the trispectrum to 
perfectly clean data. 
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Fig. 13. For the point sources: Distribution of the KS proba- 
bilities obtained by i) comparing the trispectrum estimators of 
the non-Gaussian maps to those of the Gaussian reference set 
(dashed line) , and ii) comparing the trispectrum estimators of 
the Gaussian counterparts to those of the Gaussian reference 
set (solid line). The upper panel is for the near-diagonal trispec- 
trum (r'+^) and the lower panel for the diagonal trispectrum 
estimator (T<°)). 

4.5.2. Point sources 

The histogram of probability values in figure 13 (up- 
per panel) does not show any extreme differences for 
T(+) between the Gaussian vs Gaussian and the non- 
Gaussian vs Gaussian comparisons, with only one value 
at {£, a) ~ (300, 0) having a particularly low probabil- 
ity of 10~*. The global probability for being Gaussian is 
8.8 10~*. This changes dramatically if we look at T(°) (Fig. 
13, lower panel). The diagonal estimator detects a very 
clear non-Gaussian signal, and the meta- statistics returns 
a probability of 0. It is quite surprising that the two es- 
timators show such a huge difference, although they are 
indeed completely independent, since T^+) does not have 
any Gaussian contribution, as opposed to T^^\ 

Looking at the distribution of the low-probability 
points in {i,a) space (Fig. 14 upper and lower panels), 
we find that T^^^ is localised around £ w 7000 and seems 
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Fig. 14. For the trispectrum of the point sources; Distribution, 
in {£,a) space, of the KS probabilities lower than 10^"^. 
Asterisks are for the comparison between Gaussian counter- 
parts and Gaussian reference set. Diamonds stand for the non- 
Gaussian versus Gaussian reference set comparison. The upper 
panel is for the near-diagonal trispectrum. Most diamonds are 
found at very small a in this case. The lower panel is for the 
diagonal trispectrum estimator. 



to prefer either very small or very large values of a, cor- 
responding to oblong configurations. The weak signal of 
r'"''-' seems to be mainly due to points with low £ and very 
small a (diamonds in Fig. 14 upper panel). 

4.5.3. Filaments 

For the case of the non-Gaussian filament maps, the KS 
test gives probabilities that the non-Gaussian signal is 
issued from the same process as the Gaussian reference 
set as low as a few 10~^°. Figure 15 shows the differ- 
ence in the distributions of the KS probabilities for di- 
agonal and near-diagonal trispectrum. Furthermore, the 
meta-probability obtained by applying the KS test on these 
probability distributions returns 0. In this case and simi- 
larly to the bispectrum, the non-Gaussian signatures are 
consequently undoubtedly detected. However, again the 
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Fig. 15. For the trispectrum of the filaments: Distribution 
of KS probabilities. The solid line is for the comparison of 
Gaussian counterparts vs Gaussian reference set for the near- 
diagonal estimator (T^^^). The dashed line is for the compar- 
ison of the non-Gaussian maps to reference set for the near- 
diagonal trispectrum. The dotted line is for the comparison of 
Gaussian counterparts vs reference set for the diagonal trispec- 
trum (r''^'). The dot-dashed line is for the non-Gaussian maps 
vs reference set for the diagonal estimator. 

diagonal estimator seems more sensitive than the near- 
diagonal trispectrum. As shown in Fig. 15, there are many 
more (twice as many) low probabilities (high confidence 
detection) in the diagonal case than in the near-diagonal. 

We try to visualise the morphological information con- 
tained in the trispectrum coefficients by displaying in (^, a) 
space the KS probabilities (Figs. 16 and 17, upper panels). 
Those associated with the detection of the non-Gaussian 
signatures are in diamonds, and those associated with sta- 
tistical fluctuations in the Gaussian realisations are in as- 
terisks. We plot in figures 16 and 17 the distribution of 
a and £ in bins of KS probabilities for the trispectrum 
(T*^+)) and (r*^"'). Both trispectra prefer small values of 
a (or the diagonal a ^ 2^ which correspond to similar 
configurations), which are associated with strongly elon- 
gated structures. The diagonal trispectrum exhibits sev- 
eral other low- probability "blobs" at intermediate a val- 
ues, but the corresponding histogram (the middle panel 
of Fig. 17) shows that these are somewhat less important. 
On the I axis, T^''' seems non-Gaussian mainly for values 
between 4000 and 6000, while r(+) rather detects a signal 
around £ = 2000 - 4000. All in aU, the two trispectra to- 
gether show a rich structure with clear concentrations of 
non-Gaussian signatures, not just random scatter. 

5. Conclusions 

In the present study, we investigate two of the major 
families of non-Gaussian estimators, namely Fourier-space 
based methods (the bi- and trispectrum) and wavelet- 
space based methods (the skewness and excess kurtosis 
of the wavelet coefficients). They are applied to three 
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Fig. 16. For the near-diagonal trispectrum estimator of the fil- 
aments: The KS probabilities lower than 10"'' in the {l,a) space 
(upper panel). The asterisks are for the comparison Gaussian 
counterparts versus reference set, and the diamonds for the 
comparison non-Gaussian maps versus reference set. The mid- 
dle panel shows the distribution of a for the above selected KS 
Drobabilities. The lower Danel shows the distribution of £ for 
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Fig. 17. For the diagonal trispectrum estimator of the fila- 
ments: The KS probabilities lower than 10^'^ in the (£,a) space 
(upper panel). The asterisks are for the comparison Gaussian 
counterparts versus reference set, and the diamonds for the 
comparison non-Gaussian maps versus reference set. The mid- 
dle panel shows the distribution of a for the above selected KS 
Drobabilities. The lower oanel shows the distribution of £ for 



quite different data sets chosen to represent a rather com- 
plete sample of possible non-Gaussian signatures (namely: 
point sources, filaments, non- linearly coupled signals). 
Additionally, the methods were applied to two sets of 
Gaussian realisations (called reference set and counter- 
part set) with the same power spectra as the non-Gaussian 
maps. Wc then use the Kolmogorov-Smirnov test to quan- 
tify the level of detection of the non-Gaussian signatures 
by comparing the distributions of the estimator values for 
the non-Gaussian maps versus the Gaussian reference set, 
and for the Gaussian counterpart set versus the reference 
set. The first comparison returns directly an estimate of 
the non-Gaussianity detected by a method, while the sec- 
ond case serves to illustrate the statistical fiuctuations 
within the Gaussian signal itself. Furthermore, wc have 
checked that other statistical tests such as the Kuiper and 
Anderson-Darling tests give the same results. 

We find that the filaments represent a highly non- 
Gaussian signal. This statistical character is undoubtedly 
detected by both the Fourier and wavelet based methods 
in the three and four point estimators. For the point source 
maps, the non-Gaussian signatures are very significantly 
detected with the excess kurtosis of the wavelet coefficients 
and the diagonal estimator of the trispectrum, but less sig- 
nificantly with the bispectrum and very marginally with 
the skewncss of the wavelet coefficients. This is consistent 
with the moments of the maps, where the non-Gaussianity 
shows up in the excess kurtosis as well. It is also expected 
from the physical nature of the signal: The maps contain 
the same number of negative and positive sources, on av- 
erage, which suppresses the skewncss. The type maps 
show just the inverse behaviour, in that it is now the three 
point tests (skewness of the pixel distribution and of the 
wavelet coefficients as well as the bispectrum) which de- 
tect the non-Gaussian character best. As mentioned ear- 
lier this is related to the signal itself which is constructed 
by adding a positive contribution. 

When assessing the relative level of detection, we no- 
tice that the wavelet analysis finds the non-Gaussian sig- 
natures in the maps with a non-linear coupling coeffi- 
cient of 1% at a very high confidence but only at the first 
decomposition scale (3 arc minute scale), while the bispec- 
trum detection is somewhere between 3 and 4cr. Wavelets 
are therefore better at placing constraints on /nl j a result 
which confirms the conclusion of a recent analysis of the 
COBE-DMR maps with wavelets by Cayon et al.(2002). 
The point source maps are clearly recognised as non- 
Gaussian by all methods. Finally we note that the diagonal 
trispectrum estimator is surprisingly much more sensitive 
than the nearly diagonal estimator, even though the lat- 
ter does not contain a Gaussian contribution. This shows 
that care must be taken if not all components of a given 
estimator are computed, as the signal may be extremely 
localised. A forthcoming study will address this question 
in more detail. 

Comparing the methods on a more fundamental level, 
we find as an additional advantage of the wavelet-based 
approach that it allows us to associate the non-Gaussian 
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signatures with the features that have caused them in the 
map, at all the scales where they exist. Wavelet decom- 
position thus permits to take benefit of any scale-scale 
correlations that might exist in the non-Gaussian signal. 
This can be done directly in the wavelet space or by select- 
ing the coefficients and reconstructing the signal with the 
inverse wavelet transform. These aspects have not been 
investigated in the present study. 

However, the wavelet based estimators are lacking one 
important property of the Fourier approach. The bi- and 
trispectrum can be analytically predicted without resort- 
ing to Monte Carlo based methods. They can thus be di- 
rectly related to physical phenomena, which explains their 
interest for cosmology. Additionally, the bi- and trispec- 
trum allow us to probe a very large number of geometrical 
configurations for the triplets and quadruplets as com- 
pared to the limited number of decomposition scales in- 
vestigated by the wavelets. 

The present comparison can be extended in terms of 
general strategy proposed to the groups analysing large 
CMB data sets. If we try to condense it into a few recom- 
mendations, we end up with the following steps: 

1. Look at the higher order moments in pixel space and 
the CPF first. The cumulative probability function is 
very easy and fast to calculate. Nonetheless, it seems 
to have a similar power as the bi- and trispectrum 
where just the detection of non-Gaussian signatures is 
concerned. It therefore provides a quick check if one 
can expect to find a strong signal with the more so- 
phisticated analysis techniques. We recommend that 
the maps are deconvolved with their power spectrum 
before testing the CPF. 

2. Apply wavelet based methods second. Wavelets have the 
highest detection power of the methods presented here, 
and are only moderately more expensive to apply than 
the CPF test. If indeed a signal is detected, they can 
reconstruct its origin in pixel space. 

3. Use Fourier methods for a detailed analysis. Together, 
the bi- and trispectrum furnish a large number of co- 
efficients associated with different geometrical config- 
urations. As the number of theoretical predictions for 
the bi- and trispectrum from various sources of non- 
Gaussian signatures is increasing, this allows a pre- 
cise characterisation of any strong and detected non- 
Gaussianity in a map, and the use of accurate matched 
filters to determine its possible origin. 

Clearly, combining both wavelet and Fourier based analy- 
ses seems the best strategy for studying the non-Gaussian 
signals. Together, they allow for very significant detections 
of the non-Gaussian signatures and a large amount of spa- 
tial information. In addition, they offer a theoretical basis 
to relate the physical processes to the measured quantities. 
Therefore, these methods offer currently the best chances 
of success when trying to identify the origin of the various 
sources of non-Gaussian signatures that are expected to 
contribute to the CMB signal. 
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